using System;
using System.Collections;
using System.Collections.Generic;
using NUnit.Framework;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.LinearAlgebra.IO;
using dnAnalytics.Exceptions;

namespace dnAnalytics.UnitTests.LinearAlgebra.Decomposition
{
    public abstract class AbstractQRTest
    {
        public abstract void GetDecompositionClass(String file, out Matrix matrix, out QR qr);
        private QR squareQr;
        private QR singularQr;
        private QR wideQr;
        private QR tallQr;

        private Matrix squareMatrix;
        private Matrix singularMatrix;
        private Matrix wideMatrix;
        private Matrix tallMatrix;

        public abstract Vector GetVector(int order);
        public abstract Matrix GetMatrix(int order);
        public abstract Matrix GetMatrix(int rows, int columns);

        [TestFixtureSetUp]
        public void SetUp()
        {
            MatrixMarketReader mmr = new MatrixMarketReader();
            GetDecompositionClass("../../TestMatrices/random_real_general_array_100.matrix_market", out squareMatrix, out squareQr);
            GetDecompositionClass("../../TestMatrices/gear_integer_general_coordinate_100.matrix_market", out singularMatrix, out singularQr);
            GetDecompositionClass("../../TestMatrices/random_real_general_array_10_20.matrix_market", out wideMatrix, out wideQr);
            GetDecompositionClass("../../TestMatrices/random_real_general_array_20_10.matrix_market", out tallMatrix, out tallQr);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void QRConstructorNull()
        {
            Matrix matrix = null;
            QR lu = new QR(matrix);
        }

        [Test]
        public void NotFullRank()
        {
            Assert.IsFalse(singularQr.IsFullRank());
        }

        [Test]
        public void FullRank()
        {
            Assert.IsTrue(squareQr.IsFullRank());
            Assert.IsTrue(wideQr.IsFullRank());
            Assert.IsTrue(tallQr.IsFullRank());
        }

        [Test]
        public void Factors()
        {
            Matrix q = squareQr.Q();
            Matrix r = squareQr.R();

            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
            q.Multiply(r, result);
            double error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }


            q = tallQr.Q();
            r = tallQr.R();

            result = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
            q.Multiply(r, result);

            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }

            q = wideQr.Q();
            r = wideQr.R();

            result = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
            q.Multiply(r, result);
            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }

            q = singularQr.Q();
            r = singularQr.R();
            result = GetMatrix(singularMatrix.Rows, singularMatrix.Columns);
            q.Multiply(r, result);
            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        public void FactorsResults()
        {
            Matrix q = GetMatrix(squareMatrix.Rows, squareMatrix.Rows);
            squareQr.Q(q);
            Matrix r = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
            squareQr.R(r);

            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Columns);
            q.Multiply(r, result);
            double error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - squareMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - squareMatrix[i, j]) / squareMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }


            q = GetMatrix(tallMatrix.Rows, tallMatrix.Rows);
            tallQr.Q(q);
            r = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
            tallQr.R(r);

            result = GetMatrix(tallMatrix.Rows, tallMatrix.Columns);
            q.Multiply(r, result);
            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(tallMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - tallMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - tallMatrix[i, j]) / tallMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }

            q = GetMatrix(wideMatrix.Rows, wideMatrix.Rows);
            wideQr.Q(q);
            r = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
            wideQr.R(r);

            result = GetMatrix(wideMatrix.Rows, wideMatrix.Columns);
            q.Multiply(r, result);
            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - wideMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - wideMatrix[i, j]) / wideMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }

            q = singularQr.Q();
            r = singularQr.R();
            result = GetMatrix(singularMatrix.Rows, singularMatrix.Columns);
            q.Multiply(r, result);
            error = 0;
            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    error = 0;
                    if (System.Math.Abs(squareMatrix[i, j]) < 1e10)
                    {
                        error = System.Math.Abs(result[i, j] - singularMatrix[i, j]);
                    }
                    else
                    {
                        error = System.Math.Abs((result[i, j] - singularMatrix[i, j]) / singularMatrix[i, j]);
                    }
                    Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void QResultNull()
        {
            Matrix result = null;
            squareQr.Q(result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void RResultNull()
        {
            Matrix result = null;
            squareQr.R(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void QResultNotConformable1()
        {
            Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
            squareQr.Q(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void QResultNotConformable2()
        {
            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Rows + 1);
            squareQr.Q(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void RResultNotConformable1()
        {
            Matrix result = GetMatrix(squareMatrix.Rows + 1, squareMatrix.Rows);
            squareQr.R(result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void RResultNotConformable2()
        {
            Matrix result = GetMatrix(squareMatrix.Rows, squareMatrix.Rows + 1);
            squareQr.R(result);
        }

        [Test]
        public void Determinant()
        {
            Assert.AreEqual(0, singularQr.Determinant(), Constants.ACCEPTABLE_ERROR);
            double expected = 5.006656015927541e248;
            double actual = squareQr.Determinant();
            double error = System.Math.Abs((actual - expected) / expected);
            Assert.IsTrue(error < Constants.ACCEPTABLE_ERROR);
        }

        [Test]
        public void SolveSquareMatrix()
        {
            Matrix identity = MatrixBuilder.CreateIdentityMatrix(squareMatrix.Rows);
            Matrix result = squareQr.Solve(identity);
            result = squareMatrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        public void SolveTallMatrix()
        {
            Matrix matrix = GetMatrix(tallMatrix.Rows, 5);
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    matrix[i, j] = (i + 1.0) / (j + 2.0);
                }
            }
            Matrix result = tallQr.Solve(matrix);
            Matrix test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(test[i, j], result[i, j], Constants.ACCEPTABLE_ERROR);
                }
            }
        }
        
        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNull()
        {
            Matrix matrix = null;
            squareQr.Solve(matrix);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformable()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows + 1);
            squareQr.Solve(matrix);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveMatrixInvalidOperation()
        {
            Matrix matrix = GetMatrix(singularMatrix.Rows);
            singularQr.Solve(matrix);
        }

        [Test]
        public void SolveSquareMatrixResult()
        {
            Matrix identity = MatrixBuilder.CreateIdentityMatrix(squareMatrix.Rows);
            Matrix result = GetMatrix(squareMatrix.Rows);
            squareQr.Solve(identity, result);
            result = squareMatrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.ACCEPTABLE_ERROR);
                    }
                }
            }
        }

        [Test]
        public void SolveTallMatrixResult()
        {
            Matrix matrix = GetMatrix(tallMatrix.Rows, 5);
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    matrix[i, j] = (i + 1.0) / (j + 2.0);
                }
            }
            Matrix result = GetMatrix(tallMatrix.Columns, 5);
            tallQr.Solve(matrix, result);
            Matrix test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(test[i, j], result[i, j], Constants.ACCEPTABLE_ERROR);
                }
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult1()
        {
            Matrix matrix = null;
            Matrix result = GetMatrix(squareMatrix.Rows);
            squareQr.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveMatrixNullResult2()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows);
            Matrix result = null;
            squareQr.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult1()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows + 1);
            Matrix result = GetMatrix(squareMatrix.Rows);
            squareQr.Solve(matrix, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveMatrixNotConformableResult2()
        {
            Matrix matrix = GetMatrix(squareMatrix.Rows);
            Matrix result = GetMatrix(squareMatrix.Rows + 1);
            squareQr.Solve(matrix, result);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveMatrixInvalidOperationResult()
        {
            Matrix matrix = GetMatrix(singularMatrix.Rows);
            Matrix result = GetMatrix(squareMatrix.Rows);
            singularQr.Solve(matrix, result);
        }

        [Test]
        public void SolveSquareVector()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = squareQr.Solve(vector);
            result = squareMatrix * result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0) / (i + 2.0), result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        public void SolveTallVector()
        {
            Vector vector = GetVector(tallMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = tallQr.Solve(vector);
            Vector test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * vector;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual(test[i], result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNull()
        {
            Vector vector = null;
            squareQr.Solve(vector);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformable()
        {
            Vector vector = GetVector(squareMatrix.Rows + 1);
            squareQr.Solve(vector);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveVectorInvalidOperation()
        {
            Vector vector = GetVector(singularMatrix.Rows);
            singularQr.Solve(vector);
        }

        [Test]
        public void SolveSquareVectorResult()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = GetVector(squareMatrix.Rows);
            squareQr.Solve(vector, result);
            result = squareMatrix * result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0) / (i + 2.0), result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        public void SolveTallVectorResult()
        {
            Vector vector = GetVector(tallMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }
            Vector result = GetVector(tallMatrix.Columns);
            tallQr.Solve(vector, result);
            Vector test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * vector;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual(test[i], result[i], Constants.ACCEPTABLE_ERROR);
            }
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult1()
        {
            Vector vector = null;
            Vector result = GetVector(squareMatrix.Rows);
            squareQr.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(ArgumentNullException))]
        public void SolveVectorNullResult2()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            Vector result = null;
            squareQr.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult1()
        {
            Vector vector = GetVector(squareMatrix.Rows + 1);
            Vector result = GetVector(squareMatrix.Rows);
            squareQr.Solve(vector, result);
        }

        [Test]
        [ExpectedException(typeof(NotConformableException))]
        public void SolveVectorNotConformableResult2()
        {
            Vector vector = GetVector(squareMatrix.Rows);
            Vector result = GetVector(squareMatrix.Rows + 1);
            squareQr.Solve(vector, result);
        }

        [ExpectedException(typeof(InvalidOperationException))]
        public void SolveVectorInvalidOperationResult()
        {
            Vector vector = GetVector(singularMatrix.Rows);
            Vector result = GetVector(squareMatrix.Rows);
            singularQr.Solve(vector, result);
        }
    }
}
